#version 300 es

// 设置精度
precision highp float;
precision highp int;
precision highp sampler3D;

// 输入参数
uniform vec2 uWindowSize;                    // 窗口大小
uniform bool uSampleCloudNoise;              // 是否采样云噪声
uniform float uT;                            // 时间
uniform float uTModded;                      // 修改后的时间
uniform float uRCloud0;                      // 云层底部半径
uniform float uRCloud1;                      // 云层顶部半径
uniform float uR1;                           // 大气层半径
uniform float uRSun;                         // 太阳半径
uniform float uRPlanet;                      // 行星半径
uniform vec3 uLightPos;                      // 光源位置
uniform float uMarchSamples;                 // 步进采样数
uniform float uSkyRayleighScatteringMul;     // 天空瑞利散射乘数
uniform float uRayleighScaleDiv;             // 瑞利散射比例除数
uniform float uMieScaleDiv;                  // 米氏散射比例除数
uniform float uSkyAbsorptionMul;             // 天空吸收乘数
uniform vec3 uSkyRayleighScattering;         // 天空瑞利散射系数
uniform vec3 uMieScattering;                 // 米氏散射系数
uniform vec3 uMieAbsorbtion;                 // 米氏吸收系数
uniform vec3 uCloudScattering;               // 云层散射系数
uniform float uCloudTexScale;                // 云纹理缩放
uniform float uErosionTexScale;              // 侵蚀纹理缩放
uniform float uErosionThreshold;             // 侵蚀阈值
uniform float uCloudDensityMul;              // 云密度乘数

// 纹理采样器
uniform sampler3D uCloudTex;                 // 云纹理
uniform sampler3D uErosionTex;               // 侵蚀纹理

// 输入和输出变量
in vec2 vTexCoord;                           // 纹理坐标
layout (location=0) out vec4 color;          // 输出颜色

// 数学常量
const float PI = 3.14159265359;
const float PI2 = PI * 2.;
const float inv4PI = 1. / (4. * PI);
const float PI_OVER_180 = PI / 180.0;
const float COT_HALF_FOV = 1. / (tan((30.) * PI_OVER_180));
const vec3 PLANET_ORIGIN = vec3(0.);

// 大气参数
const float RAYLEIGH_SCATTERING_MULTIPLIER = .03624;    // 瑞利散射乘数
const float MIE_SCATTERING_MULTIPLIER = .00692;         // 米氏散射乘数
const float MIE_ABSORBTION_MULTIPLIER = .00077;         // 米氏吸收乘数
const float OZONE_ABSORBTION_LAYER_WIDTH = 25.;         // 臭氧吸收层宽度
const float OZONE_LINEAR_TERM0 = 1. / 15.;              // 臭氧线性项0
const float OZONE_CONSTANT_TERM0 = -2. / 3.;            // 臭氧常数项0
const float OZONE_LINEAR_TERM1 = -1. / 15.;             // 臭氧线性项1
const float OZONE_CONSTANT_TERM1 = 8. / 3.;             // 臭氧常数项1
const vec3 OZONE_ABSORBTION = vec3(.000650, .001881, .000085);  // 臭氧吸收系数

/**
 * 阴影步进数据结构体
 * 用于存储阴影计算过程中的数据
 */
struct ShadowMarchData {
  vec3 opticalDepth;  // 光学深度
  vec3 nearSample;    // 近采样点
  vec3 farSample;     // 远采样点
};

/**
 * 参与介质结构体
 * 用于描述介质的光学特性
 */
struct ParticipatingMedia {
  vec3 rayleighScattering;  // 瑞利散射系数
  vec3 mieScattering;       // 米氏散射系数
  vec3 scattering;          // 总散射系数
  vec3 extinction;          // 消光系数
};

/**
 * 生成随机数
 * @param n - 随机数种子
 * @return float - 随机数
 */
float nrand(vec2 n) {
  return fract(sin(dot(n.xy, vec2(12.9898, 78.233)))* 43758.5453);
}

/**
 * 生成带时间变化的随机数
 * @param n - 随机数种子
 * @return float - 随机数
 */
float n1rand(vec2 n) {
  float t = fract(uTModded);
  float nrnd0 = nrand(n + 0.07 * t);
  return nrnd0;
}

/**
 * 计算射线与球体的交点
 * @param r0 - 射线起点
 * @param rd - 射线方向
 * @param s0 - 球心
 * @param sR - 球半径
 * @param farthest - 是否返回最远交点
 * @return float - 交点距离，-1表示无交点
 */
float raySphereIntersect(vec3 r0, vec3 rd, vec3 s0, float sR, bool farthest) {
  float a = dot(rd, rd);
  vec3 s0_r0 = r0 - s0;
  float b = 2.0 * dot(rd, s0_r0);
  float c = dot(s0_r0, s0_r0) - (sR * sR);
  float delta = b * b - 4.0*a*c;
  if (delta < 0.0 || a == 0.0) {
    return -1.0;
  }
  float sol0 = (-b - sqrt(delta)) / (2.0*a);
  float sol1 = (-b + sqrt(delta)) / (2.0*a);
  if (sol0 < 0.0 && sol1 < 0.0) {
    return -1.0;
  }
  if (sol0 < 0.0) {
    return max(0.0, sol1);
  } else if (sol1 < 0.0) {
    return max(0.0, sol0);
  }

  if (farthest) {
    return max(0.0, max(sol0, sol1));
  } else {
    return max(0.0, min(sol0, sol1));
  }
}

/**
 * 将UV坐标转换为透射率参数
 * @param viewHeight - 输出视图高度
 * @param viewZenithCosAngle - 输出天顶角余弦
 * @param uv - 输入UV坐标
 */
void uvToLutTransmittanceParams (out float viewHeight, out float viewZenithCosAngle, in vec2 uv) {
  float x_mu = uv.x;
  float x_r = uv.y;

  float H = sqrt(uR1 * uR1 - uRPlanet * uRPlanet);
  float rho = H * x_r;
  viewHeight = sqrt(rho * rho + uRPlanet * uRPlanet);

  float d_min = uR1 - viewHeight;
  float d_max = rho + H;
  float d = d_min + x_mu * (d_max - d_min);
  viewZenithCosAngle = d == 0.0 ? 1.0f : (H * H - rho * rho - d * d) / (2.0 * viewHeight * d);
  viewZenithCosAngle = clamp(viewZenithCosAngle, -1.0, 1.0);
}

/**
 * 获取天空参与介质
 * @param samplePos - 采样位置
 * @param RAYLEIGH_EXP_SCALE - 瑞利散射指数比例
 * @param MIE_EXP_SCALE - 米氏散射指数比例
 * @return ParticipatingMedia - 参与介质数据
 */
ParticipatingMedia getSkyParticipatingMedia(const vec3 samplePos, const float RAYLEIGH_EXP_SCALE, const float MIE_EXP_SCALE) {
  float sampleAltitude = length(samplePos) - uRPlanet;
  float rayleightDensity = exp(-sampleAltitude * RAYLEIGH_EXP_SCALE);
  float mieDensity = exp(-sampleAltitude * MIE_EXP_SCALE);
  float ozoneDensity = clamp(sampleAltitude < OZONE_ABSORBTION_LAYER_WIDTH ? OZONE_LINEAR_TERM0 * sampleAltitude + OZONE_CONSTANT_TERM0 : OZONE_LINEAR_TERM1 * sampleAltitude + OZONE_CONSTANT_TERM1, 0., 1.);
  vec3 mieAbsorbtion = uMieAbsorbtion * mieDensity;
  vec3 mieScattering = uMieScattering * mieDensity;
  vec3 rayleighScattering = uSkyRayleighScattering * uSkyRayleighScatteringMul * rayleightDensity;
  vec3 ozoneAbsorbtion = OZONE_ABSORBTION * ozoneDensity * uSkyAbsorptionMul;
  vec3 sigmaScattering = mieScattering + rayleighScattering;
  vec3 sigmaAbsorption = mieAbsorbtion + ozoneAbsorbtion;
  vec3 sigmaExtinction = max(vec3(.000000001), sigmaAbsorption + sigmaScattering);

  return ParticipatingMedia(rayleighScattering, mieScattering, sigmaScattering, sigmaExtinction);
}

/**
 * 步进天空层
 * @param marchOrigin - 步进起点
 * @param rayDir - 射线方向
 * @param marchDist - 步进距离
 * @param steps - 步数
 * @return ShadowMarchData - 阴影步进数据
 */
ShadowMarchData marchSkyLayer (vec3 marchOrigin, vec3 rayDir, float marchDist, float steps) {
  vec3 samplePos;
  vec3 nearSample;
  vec3 farSample;

  vec3 sigmaExtinction;
  vec3 opticalDepth = vec3(0.);

  float stepSize = marchDist / steps;
  float depth = .001;
  float RAYLEIGH_EXP_SCALE = 1. / uRayleighScaleDiv;
  float MIE_EXP_SCALE = 1. / uMieScaleDiv;

  while (depth < marchDist) {
    samplePos = marchOrigin + rayDir * depth;
    sigmaExtinction = getSkyParticipatingMedia(samplePos, RAYLEIGH_EXP_SCALE, MIE_EXP_SCALE).extinction;
    opticalDepth += sigmaExtinction * stepSize;
    depth += stepSize;
  }

  nearSample = marchOrigin;
  farSample = samplePos;

  return ShadowMarchData(opticalDepth, nearSample, farSample);
}

/**
 * 步进大气层
 * @param cameraPos - 相机位置
 * @param rayDir - 射线方向
 * @return vec4 - 渲染结果
 */
vec4 marchAtmosphere (vec3 cameraPos, vec3 rayDir) {
  float marchDist;
  float R1Dist = raySphereIntersect(cameraPos, rayDir, PLANET_ORIGIN, uR1, false);
  float RPlanetDist = raySphereIntersect(cameraPos, rayDir, PLANET_ORIGIN, uRPlanet, false);

  if (R1Dist == -1.) {
    return vec4(0.);
  }

  if (RPlanetDist == -1.) {
    marchDist = R1Dist;
  } else {
    marchDist = min(R1Dist, RPlanetDist);
  }

  ShadowMarchData skyLayer = marchSkyLayer(cameraPos, rayDir, marchDist, uMarchSamples);
  return vec4(skyLayer.opticalDepth, 1.);
}

/**
 * 主函数
 * 渲染天空阴影贴图
 */
void main () {
  // 获取UV坐标
  vec2 uv = vTexCoord;
  
  // 计算视图参数
  float viewHeight;
  float viewZenithCosAngle;
  uvToLutTransmittanceParams(viewHeight, viewZenithCosAngle, uv);
  
  // 设置相机位置和射线方向
  vec3 cameraPos = vec3(0., 0., viewHeight);
  vec3 rayDir = vec3(0., sqrt(1. - viewZenithCosAngle * viewZenithCosAngle), viewZenithCosAngle);
  
  // 计算大气层光学深度
  color = marchAtmosphere(cameraPos, rayDir);
}
